Local and non-local electron-phonon couplings in K 3 -picene and the effect of metallic 

screening. 
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We analyze the properties of electron-phonon couplings in K3-picene by exploiting a molecular 
orbital representation derived in the maximally localized Wannier function formalism. This allows 
us to go beyond the analysis done in Phys. Rev. Lett. 107, 137006 (2011), and separate not only 
the intra- and intermolecular phonon contributions but also the local and non-local electronic states 
in the electron-phonon matrix elements. Despite the molecular nature of the crystal, we find that 
the purely molecular contributions (Holstein-like couplings where the local deformation potential is 
coupled to intramolecular phonons) account for only 20% of the total electron-phonon interaction 
A. In particular, the Holstein-like contributions to A in K3-picene are four times smaller than those 
computed for an isolated neutral molecule, as they are strongly screened by the metallic bands 
of the doped crystal. Our findings invalidate the use of molecular electron-phonon calculations to 
estimate the total electron-phonon coupling in metallic picene, and possibly in other doped metallic 
molecular crystals. The major contribution (80%) to A in K3-picene comes from non-local couplings 
due to phonon modulated hoppings. We show that the crystal geometry together with the molecular 
picene structure leads to a strong ID spatial anisotropy of the non-local couplings. Finally, based 
on the parameters derived from our density functional theory calculations, we propose a lattice 
modelization of the electron-phonon couplings in K3-picene which gives 90% of ab-initio A. 

PACS numbers: 31.15.A-, 74.70.Kn, 63.20.kd 



I. INTRODUCTION 

Understanding transport properties of molecular crys- 
tals based on hydrocarbon molecules is relevant not only 
to fundamental condensed matter physics, but also for 
applications in nanoelectronics. For instance, organic 
field-effect transistors are appealing as they are flexi- 
ble, lightweight and cheap. Rubrene-based field effect 
transistors 1,2 display tunable mobilities that can be as 
large as 40 cm 2 / (V • s). More recently it was shown that 
picene field-effect transistors 3 based on liquid electrolytes 
have p-channel characteristics 4 , although with much re- 
duced mobilities with respect to rubrene. 

Transport properties of organic molecular crystals can 
also be tuned by intercalation of alkali or alkaline earth 
metals. K intercalation leads to metallic states in ph- 
thalocyanine materials 5 and in several other polycyclic 
aromatic hydrocarbons. In picene 6 , phenanthrene 7,8 , 
coronene 9 , and in l,2:8,9-dibenzopentacene 10 , intercala- 
tion stabilizes a superconducting state with critical tem- 
peratures (T c ) up to 33 K. A detailed understanding of 
transport phenomena in such systems is then relevant 
also for the realm of fundamental research. 

An important source of intrinsic scattering in aromatic 
molecular crystals is provided by the electron-phonon 
coupling. In these systems there is an interplay between 
intramolecular local interactions and intermolecular non- 
local interactions. Determining the mutual role of lo- 
cal and non-local interactions is hardly doable without a 
proper theoretical approach. Molecular crystals can in- 
deed behave very differently depending on the details of 
the molecules composing the crystal and on their arrange- 



ment. In alkali doped fullerenes 11 , superconductivity is 
supposed to be mostly due to intramolecular phonons. In 
this case, if the electronic states coupled to the phonons 
are molecular and the metallic screening is weak, then 
the problem can be tackled at a molecular level by the 
calculation of electron-phonon interaction 12 ' 13 on an iso- 
lated ionized molecule. 

The situation is more complicated in the field of hy- 
drocarbon molecular crystals. In the case of K 3 -picene, 
molecular calculations 14 using the B3LYP functional give 
a large electron-phonon coupling, that can almost alone 
explain T c . However, the generalization of this ap- 
proach to other hydrocarbon molecular crystals predicts 
a decrease of the critical temperature with the increase 
of the molecular size, in disagreement with experimen- 
tal data 6-8 ' 10 . Indeed in experiment the largest T c 
is for the crystal composed by the largest molecules. 
Subedi et al. 15 performed a density functional theory 
calculations (DFT) in which the crystal structure of pris- 
tine picene was adopted and K-doping was treated in 
a rigid doping approach. The screening of the self- 
consistent potential was assumed to be that of insulating 
picene. A very large electron-phonon coupling was found 
mostly due to intramolecular phonons, in agreement with 
Ref. 14. 

In our previous work 16 we performed DFT calcula- 
tions relying on less approximations then in Ref. 15. The 
theoretically devised crystal structure of K3-picene was 
considered 16,17 and K-atoms were explicitly included in 
the calculation. Furthermore we included the metallic 
screening of crystalline K3-picene in the self-consistent 
potential. By projecting the phonon polarizations into 
intramolecular and intermolecular vibrations, we found 
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that K3-picene has a strong electron-phonon coupling 
(A = 0.73) that is partially due to the coupling to inter- 
molecular and intercalant phonons (40%) and partially 
to the coupling to intramolecular phonons (60%), in dis- 
agreement with Refs. 14,15. In our present work, we go 
beyond what we have done in Ref. 16. Instead of ana- 
lyzing the "locality" of the electron-phonon coupling in 
terms of phonon projections only, we study it also by 
means of electronic projections onto a molecular basis, 
which allows one to distinguish between the on-site elec- 
tronic Hamiltonian and the hopping parts, both modu- 
lated by the coupling with phonons. This approach leads 
to a stricter distinction between purely molecular and 
crystal contributions, and yields a further reduction of 
the purely molecular component, estimated to be about 
20% of the total A. 

The three approaches illustrated above, namely molec- 
ular calculations, rigid doping of the crystal, and explicit 
treatment of the dopants, rely on different approxima- 
tions that could explain the discrepancies. An important 
one is the treatment of the electronic screening and its 
effects on the electron-phonon interaction. In molecular 
calculations 14 and in Ref. 15, metallic screening is ne- 
glected. Analogy with alkali doped fullerenes points out, 
however, that this assumption is not necessary fulfilled. 
In K3C60, it has been suggested that metallic screening 
strongly affects the electron-phonon coupling 18-20 . For 
example, Ai s modes causing a shift without splitting of 
the ti u Geo molecular levels, are supposed to be screened 
by the charge transfer from up-shifted to down-shifted 
levels 20 , i. e. by the metallic screening in the solid. 
In K 3 -picene the situation could be similar. However, 
the relative contribution of intramolecular, intermolec- 
ular and intercalant interactions remains unclear, and 
largely unexplored. In this work we carry out a detailed 
and quantitative analysis of the total electron-phonon 
coupling A in K 3 -picene, by addressing these issues. 

The paper is organized as follows. In Sec. II we pro- 
vide the general definition of local and non-local electron- 
phonon couplings for a molecular crystal. In Sec. Ill we 
describe the geometry of K3-picene, and we show the 
dominant hoppings of the corresponding tight-binding 
Hamiltonian in the Wannier basis. In Sec. IV electron- 
phonon calculations are carried out by discriminating be- 
tween local and non-local couplings. We find that the 
purely local contributions account for only 20% of the full 
A, while the remaining part comes from non-local sources. 
Sec V analyzes the screening acting on the local electron- 
phonon terms by a direct comparison between the crystal 
and the isolated (unscreened) molecule. We show that 
the effect of the metallic screening provided by the crys- 
tal environment to the deformation potential is sizable 
with a strong reduction of the local electron-phonon cou- 
pling with respect to the corresponding strength found in 
the neutral isolated molecule. In Sec. VI we look for the 
most important non-local terms contributing to the total 
A and we build a model Hamiltonian with few non-local 
electron-phonon couplings added to the local part which 



gives 90% of the total A. The conclusions are in Sec. VII. 

II. DEFINITION OF LOCAL AND NON-LOCAL 
ELECTRON-PHONON COUPLINGS IN A 
MOLECULAR CRYSTAL 

We suppose that the band structure of a molecular 
crystal is described by the electronic tight-binding Hamil- 
tonian H c \ written in a basis set built out of molecular 
orbitals \i,m) — c| TO |0), where i is the index of the molec- 
ular site having its center of mass located at the equilib- 
rium position Hi, and m is the orbital index, with c\ m 
and Cj n satisfying canonical anticommutation relations. 
For simplicity, we assume here that there is only one 
molecule per unit cell, so the vectors Ri define also the 
Bravais lattice. H c \ reads then as 

H c\ = - ^2 ^2 t mn (Rj - Rj)c im Cj„, (1) 

ij mn 

where we omitted the spin index by implicitly assuming 
that the spin up and spin down components are equiv- 
alent, namely there is no spin symmetry breaking. The 
hopping matrix is defined as: 

-t mn (TL j --R i ) = (i,m\H\j,n), (2) 

where we exploit the lattice translational invariance. 

In second quantization a phonon displacement u s of 
atom s with mass M s relative to the i-th molecule is: 

(3) 

where N q is the number of phonon momentum points 
describing the system, Wqy is the phonon dispersion of 
mode v at a given momentum q, e q „ is the 3-dimensional 
s-atomic component of the phonon eigenvector e q „, and 
t s is the position of the atom s in the unit cell. The oper- 
ators & q „ and bqv satisfy canonical bosonic commutation 
relations. 

The harmonic phonon Hamiltonian H p h on reads as: 




The complete Hamiltonian for the electron-phonon (el- 
phon) problem includes electron-phonon coupling terms, 
and is written as 

tt tt 1 tt 1 rrlocal 1 n-non local /r\ 

£1 — n e \ -t- -TJphon i- ^d-phon "+" -"cl-phon \°) 

where the local electron-phonon coupling is 

^ i mn 

(6) 
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while the non-local coupling #ei-phon al 1S 

r E E E - R«) e tq ' R *cL Cj „6 q , + h. c .) . 



A, 



The phase e jq Ri makes the total momentum conserved 
in the electron-phonon scattering terms. The electron- 
phonon coupling strength projected on the molecular or- 
bitals is defined as 

(/^(Rj-Ri) = y;<», m\-^-\j, n)-e'„/ V 2M « w q«" (8) 



q,s 



where z,j are indexes of molecular sites, u qs is the 
Fourier transform of the phonon displacement u s (Rj), 
and 5v/Suq S is the (screened) deformation potential. 

From Eq. 6, it is apparent that the local coupling is a 
Holstein-type interaction which couples the phonons with 
on-site molecular electronic terms, while in Eq. 7 the non- 
local couplings modulate the hoppings t TOn (Rj — Rj) in 
H c \ via the bosonic fields &_ q „ and bq V . The local and 
non-local coupling strengths are proportional to the de- 
formation potential expressed in the molecular orbital 
basis, centered on cither the same site or two different 
molecules, respectively. By translational invariance, the 
strength depends only on the vector Rj — R^. It is worth 
pointing out that in this context the definition of "lo- 
cal" and "non-local" couplings is purely electronic. In 
our previous work 16 , we distinguished between the "inter- 
molecular" and "intramolecular" contributions based on 
the phonon projections. The intramolecular phonons are 
those having e q „ projected on the single molecule many- 
fold, while the intermolecular phonons are those having 
e q „ spanned by the rigid molecular rototranslations to- 
gether with all intercalant displacements. Therefore, one 
can expand the bosonic fields 6 qj , into 6 q * or + & q * ra , be- 
ing the sum of intermolecular and intramolecular pro- 
jections a resolution of the identity. We thus note that 
cross-contributions like intramolecular phonons in non- 
local couplings or intermolecular phonons in local cou- 
plings, are possible. Projection of both the electronic 
and phononic parts guarantees the isolation of the single 
molecule contribution. In this work, we are going to use 
the words "local" and "non-local couplings" to mean the 
electronic molecular basis set projections as in Eqs. 6 and 
7, respectively, while we keep the notation of Ref. 16 by 
using "intramolecular" and "intermolecular phonons" to 
refer to the phonon projections. 



III. GEOMETRY AND BAND STRUCTURE OF 
K3PICENE 

The molecules in the K 3 -picene crystal are arranged to 
satisfy the P 2 symmetry group. The unit cell contains 
two molecules and is monoclinic with axes a = 8.707A, 
b = 5.912A, c = 12.97A, a = 90°, P = 92,77°,7 = 



90°. The unit cell parameters have been taken from the 
experiment 6 , while the internal coordinates have been op- 
timized after a full geometry relaxation performed in a 
DFT framework within the local density approximation 
(LDA) (for more details see the supplementary materials 
section of Ref. 16). The final structure is drawn in Fig. 1, 
where we plot the orthogonal projections of the unit cell 
repeated twice in each crystallographic direction. From 
Figs. 1(b) and 1(c), one can clearly see the molecular 
stacking along the c-axis, while in Fig. 1(a) (the ab pro- 
jection) the molecular herringbone arrangement of each 
layer is visible. The intercalant occupies the interstitial 
space and tunes the intermolecular angles by steric effect. 

In order to understand the interplay between the K 3 - 
picene geometry and its band structure, we are going to 
derive a tight-binding model constructed on a Wannier 
function basis. The maximally localized Wannier repre- 
sentation of the DFT orbitals is useful not only to im- 
plement an interpolation scheme for computing the band 
structure and the electron-phonon matrix elements, but 
also to have a physical insight on the system. For ex- 
ample, the formation of the chemical bond in a solid can 
be visualized by means of the Wannier representation of 
the molecular orbitals (MOs). In a molecular crystal, as 
the picene, the Wannier representation is even more nat- 
ural, as it builds on the local nature of molecular sites, 
where the MOs are strongly localized. The spatial lo- 
cal representation given by the Wannier transformation 
helps in modeling the electronic structure of the K-doped 
picene, and understanding the mechanism which sets the 
superconductivity. 

By following Ref. 22, the maximally localized Wannier 
functions (MLWF) are defined as 



1 



E 



E^ 



MLWF 



(k)V>mk(r) 



-ik-R 



( 9 ) 

where the sum ^ k is over a iV w -point grid in the Brillouin 
zone (BZ) 23 , R is a Bravais lattice vector, V'mk(r) are the 
Bloch eigenstates of the rn-th band, and XJ MLWF (k) is 
a unitary matrix (for composite bands) , defined to mini- 
mize the total spread of the wave function 



X ^ [(w n o\r 2 \w n0 ) - \(w n0 \r\w n0 )\ 2 ] 



(10) 



Note that in this case there are two molecules per unit 
cell, and so the Bravais vectors R are not the centers of 
each molecule, at variance with the simplest case taken 
into account in Sec. II. In the K 3 -picene, the MLWFs 
have been determined for the bands derived from the 
lowest unoccupied molecular orbital (LUMO), LUMO+1, 
and LUMO+2 of the neutral picene molecule. Those 
bands form a quasi-composite group, as the LUMO is 
well separated from the highest occupied molecular or- 
bital (HOMO). The HOMO-LUMO gap in the pristine 
picene is 3 eV large 24 , and only the LUMO+2 is weakly 
entangled with the upper bands. Therefore, a prelimi- 
nary disentanglement procedure has been performed 25 , 
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FIG. 1: (color online) Orthogonal projections of the K3-picene unit cell, repeated twice along each crystallographic direction. 
The unit cell sides are drawn in red. Carbon atoms are in yellow, Hydrogen is in blue, and Potassium in green. Panel (a): ab 
projection with the c-axis pointing outwards; panel (b): ac projection with the fo-axis pointing inwards; panel (c): be projection 
with the a-axis pointing inwards 21 . 



before JJ MLWF (k) could be obtained. Thus, in our case 
TJ MLWF (k) is a 6 x 6 matrix (3 bands per molecule, 2 
molecules per unit cell) , and a tight-binding Hamiltonian 
can be denned in the rotated MLWF basis, according to 
the matrix elements 

H nrn (R) — (w n o\H\w m n) , (11) 

where H is the one-body LDA Hamiltonian. 

In molecular crystals, the MLWF is not necessarily the 
best basis to work with. The most "physical" basis is 
the one which diagonalizes the local part (H nm (0) and 
(n, m) running on the same molecule) of the Hamilto- 
nian in Eq. 11. Indeed, the local part of H represents 
the molecule in the crystal, and its eigenvectors w n and 
eigenvalues e™° 1 are respectively the MOs and molec- 
ular levels in the crystal environment. From here on, 
we define w ra R,(r) to be the "molecular" MLWFs, where 
U(k) = U MLWF (k) x U mo i, with U mo i being the unitary 
transformation which diagonalizes the local problem in 
the MLWF basis. The molecular MLWFs \w nR ) obtained 
in the rigorous Wannier function formalism play the role 
of the molecular orbitals \i,m) generically introduced in 
Sec. II. 

In order to see how the crystal environment affects 
the local MOs, in Fig. 2 we plotted the local molecular 
MLWF functions w n o (r) of K 3 -picene together with the 
MOs of the isolated neutral molecule (w^ OL (r)). One 
can see that w n o (r) in the doped crystal are a good rep- 
resentation of the orbitals in the isolated neutral picene. 
Indeed, the LUMO and LUMO+1 are in close agreement. 



The LUMO+2 differs only slightly, as in the crystal it is 
more "delocalized" , something expected as in the molecu- 
lar calculations its energy level is close to the free particle 
continuum, and therefore it is more affected by the en- 
vironment. The overall agreement allows one to make a 
one-to-one correspondence between the molecular prop- 
erties and the crystal local on-site properties expressed 
in the molecular MLWF. 

We now analyze the hopping terms in Eq. 11. They 
show a clear hierarchy in magnitude depending on their 
spatial direction. The largest are the nearest neighbors 
(NN) hoppings which connect the molecules within the 
herringbone layer. In the herringbone structure, each 
molecule is linked to its four nearest neighbors in two dif- 
ferent ways, by the proximity of either a two-ring molec- 
ular side, or a three-ring side (see Fig. 3(a)). We found 
that there is a large asymmetry between the NN hoppings 
connecting two molecules via a three-ring molecular side, 
dubbed "ID NN" in the text, and the ones whose con- 
nection is bridged by a two-ring side, dubbed "2D NN" . 
The ID NN terms, sized up to 0.09 eV, are almost twice 
larger than the 2D NN hoppings, which reach 0.05 eV at 
most. This is due to the internal degrees of freedom of the 
single crystal site, as the picene has an aromatic 5-ring 
structure. If the molecule were symmetric, the ID and 
2D NN hopping terms would be equal. The consequences 
of this internal asymmetry will be studied later in both 
the band structure and the electron-phonon couplings. 

Not only the NN but also the next-nearest neigh- 
bors (NNN) hoppings are not symmetric. Indeed, the 
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molecular LUMO 



molecular LUMO+1 



molecular LUMO+2 



crystal LUMO^^ crystal LUMO+1 crystal LUMO+2 

FIG. 2: (color online) Surface plot of the orbitals of an isolated neutral picene molecule (upper row) and the molecular MLWF 
functions (ui n o(r)) in the K3-picene crystal (lower row). The surface is defined by the set of points which satisfies the condition 
|*l/n(r)| = 0.05, where f n (r) = Re[</> n (r)]Sign[</> n (r)], with 4> n {v) the MO whose phase has been fixed and normalized such that 
</>n(r max ) = 1, being r max the location of the maximum of its modulus 21 . 
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FIG. 3: (color online) Intermolecular hoppings in the herringbone layer. Panel (a): The two possible nearest-neighbors (NN) 
intermolecular hoppings are represented, with the strongest in orange (ID NN) mediated by a 3-ring molecular side, while the 
weakest in blue (2D NN) is bridged by a 2-ring side. Panel (b): All possible NN and NNN hoppings in the plane are drawn. 
The ID NNN are in red, while the 2D NNN are in dark blue. Panel (c): two ladder chains formed by selecting only the ID NN 
and ID NNN hoppings, the strongest ones among the bidimensional hoppings. 



NNN terms pointing along the b crystallographic axis 
(named "ID NNN" in the text) are more than twice 
larger than the ones pointing along the a crystallographic 
axis (dubbed "2D NNN"), which do not go beyond 0.02 
cV. This can be easily explained by noting that the b 
axis is shorter than the a axis, thus in the b direction the 
molecules are more closely packed, with an increase of the 
transfer integrals and so of the hoppings. See Fig. 3(b) 
for the graphical representation of all the NN and NNN 
hoppings in the herringbone plane. 

It turns out that the ID NNN and the 2D NN terms are 
of the same magnitude (sa 0.05 eV). The combination of 
ID NN and ID NNN hoppings only, creates ladder chains 
spanning the b axis (see Fig. 3(c)), while in a four hop- 
ping model (with the addition of the 2D NN and the 2D 
NNN terms), their combination spans the full 2D space. 
The 2-ring versus 3-ring asymmetry clearly favors a "ne- 
matic" one dimensional electronic structure with respect 
to the full bidimensional layer. Therefore we define a 



"ID model" comprising of the 0D (on-site), ID NN and 
ID NNN terms, and a "2D model" which includes all 
terms of the "ID model" plus the 2D NN and 2D NNN 
hoppings. 

To understand the impact of this hierarchy on the 
band structure we take the hoppings of the tight-binding 
Hamiltonian written in the MLWF basis, and we are 
going to selectively switch them on and off. The full 
band structure is plotted in Fig. 4(d) for the LUMO, 
LUMO+1, and LUMO+2 states, which yield 2 bands 
each. By keeping only the local on-site terms, we ob- 
tain the molecular levels e„ in the crystal, which of 
course are dispersionless (Fig. 4(a)). By switching on 
the ID NN and the ID NNN hoppings along the molecu- 
lar "wire" , one gets the band structure of the ID model 
in Fig. 4(b). The dispersion develops only along the 
b-axis, but it gives the main contribution to the full 3D 
bandwidth, while the double degeneracy along the CY 
path is due to the P-j 2 symmetry. The full band struc- 
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ture can be roughly modeled by the ID model, except 
that the Fermi surface is poorly reproduced. For the 2D 
model (see Fig. 1(b)), a band structure closer to the 3D 
one is obtained, with the LUMO bands almost perfectly 
reproduced, and the flatness of the BD and CY paths 
due to the decoupling in the layer stacking. 




. (a) 3D without LUMO+2 

- 


. (b) 3D 















r B D Z CY r B D Z CY r 

FIG. 5: The full band structure (right) and the one obtained 
by taking off the LUMO+2 states (left) from the tight-binding 
Hamiltonian obtained in the molecular MLWF basis for the 
K3-picene. The definition of the high-symmetry points in the 
k-path is reported in the caption of Fig. 4. The zero of the 
energy axis is the Fermi level. 



Technical details for the band structure 
calculations 



r B D z 



C Y 



r B D z 



C Y 



FIG. 4: Band structure of a tight-binding Hamiltonian for 
the K3-picene derived in a MLWF basis including LUMO, 
LUMO+1 and LUMO+2 states, plotted along the TBDZCYr 
k-space path. In the reciprocal crystal units, the special 
points arc T = (0,0,0), B= (£,0,0), D= (i,0, ±), Z= 
(0,0, |), C= (0, |,i), Y= (0, |,0). In the left-upper panel 
only the "on-site" hoppings have been retained, while in the 
right-lower panel the band structure has been obtained with 
the full tight-binding model. The upper-right (the lower-left) 
panel is the result of a tight-binding model with only ID 
(1D+2D) nearest neighbor and next-nearest neighbor hop- 
pings. The zero of the energy axis is the Fermi level. 



In the spirit of downfolding the full electronic structure 
to a low-energy lattice model, one interesting question 
is whether a two-orbital model is enough to reproduce 
the low-energy physics. To this aim, we suppressed the 
LUMO+2 orbital from the tight-binding model. The re- 
sult is shown in Fig. 5. As one can see, the LUMO+1 
bands are strongly deformed, and the Fermi surface is 
strongly modified. To have a correct description of the 
low-energy physics of the crystal, one needs also to in- 
clude the LUMO+2 molecular orbital. Therefore, a cor- 
rect modclization of the system comprises 3 orbitals, up 
to the LUMO+2, being the hybridization between the 
LUMO+1 and LUMO+2 very strong. 



The LDA-DFT calculations have been performed with 
the Quantum-espresso 26 code. K, C, and H atoms are 
described by ultrasoft pseudopotentials. The plane-wave 
(PW) cutoff is 60 Ry for the wave- function, and 600 Ry 
for the charge. A 4 x 4 x 4 electron-momentum grid and a 
Methfessel-Paxton smearing of 0.015 Ry are used in the 
electronic integration. 

Wannierization has been performed with the 
Wannier90 27 program on a N w =4x4x4 electron- 
momentum mesh, by including the LUMO, LUMO+1, 
and LUMO+2 states. Both the long-range hoppings and 
inclusion of the first three LUMO's are needed to get 
localized orbitals and Wannierized bands in a very good 
agreement with those computed in the PW basis set in a 
window of ± 0.3 eV around the Fermi level. Indeed, the 
maximum discrepancy between the ab-initio bands and 
the Wannierized ones is only 5 meV for the LUMO and 
LUMO+1 states, while it is larger (0.05 eV at most) for 
the LUMO+2 band, that is however higher in energy. 

The spreads fl n = (w nQ \r 2 \w n0 ) - |(w n o|r|*no)| 2 of 
the n-th MLWF w nR (r) are 11.8A, 13.2A, and 21.1 Afor 
n — 1, 2, and 3, respectively. The n — 3 MLWF orbital 
is more spread out, and leads to a LUMO+2 molecu- 
lar orbital more sensitive to the crystal environment, as 
highlighted by Fig. 2. 
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IV. LOCAL AND NON-LOCAL 
ELECTRON-PHONON COUPLINGS IN 
K3PICENE 



The 

iVo E 



total electron-phonon coupling is A = 
ql/ A q ^, where v is the phonon mode and q is 
its momentum. The phonon resolved coupling reads: 



W 2 JV(0) A k 



k+qm I 



qy \ j k „ jTO 

x (/kn — /k+q,m) <$(ek+q,m — £kn — W q „), (12) 

that couples the occupied state |k, n) (the ket refers 
to the periodic part of the Bloch function) of momen- 
tum k and band n with the empty state |k + q, m) 
separated by the phonon energy w qi/ . N(0) is the 
electron DOS per spin per cell at the Fermi level. 
The electron-phonon matrix elements are g£ n k+qm = 
E s e q ,-d^ n (k+q, k)/ v /2M^ q 7, where d^„(k+q, k) - 
(k + q, m\dvscF / Suq S \k, n), with Sv SC F/Suq S the peri- 
odic part of the DFT screened deformation potential. In 
Eq. 12, /k„ are Fermi functions depending on the temper- 
ature T, and the expression for A q „ has to be evaluated by 
a T — > extrapolation. In the "adiabatic" limit, namely 
for ujph <C Ae, where Ae is the bandwidth and uj p h is the 
characteristic phonon frequency, the expression for X qu 
in Eq. 12 reduces to the one proposed by Allen 28 , and 
generally used in previous electron-phonon estimates: 



J lfkn,k+qn 



\AD ± 

X <5(ek,n)<5(ek+q,m)- 



(13) 



We are going to dub X AD in the Equation above as "adi- 
abatic", while A in Eq. 12 as "non-adiabatic" . 

By exploiting the definition of Wannier functions 
in Eq. 9, the deformation potential matrix elements 
d^ nn (k + q, k) can be written in terms of the molecular 
MWLF basis as 



.ik R 



C( k + q> k ) = ErE„ 

t/ roro ,(k + q)d^ n ,(R)£0(k), 



(14) 



where the deformation potential in the MLWF local rep- 
resentation is 



d mn( R ) = (^mo| ^ SCF |w»R)- 



5u. 



qs 



Eq. 15 is the analogous of Eq. 11 but for the electron- 
phonon coupling elements, when only the localization of 
the wave function is used. Therefore, the same analysis 
carried out in Sec. Ill can be done here, with the distinc- 
tion between the "local" (with R = and (to, n) orbitals 
on a single molecule) and "non-local" (with R ^ 0, or 
R = with (m,n) orbitals centered on two different 
molecules of the unit cell) matrix elements. Thus, it 
is the Wannier function formalism which allows one to 



make the bridge from the plane wave representation to 
the molecular orbital description of the electron-phonon 
problem introduced in Sec. II, with the distinction be- 
tween local and non-local couplings. 

As already pointed out in Sec. II, an analogous but in- 
dependent definition of local and non-local contributions 
can be done not only for the electronic states, but also for 
the phonon modes. To project the phonon vibrations we 
use the same strategy as reported in Ref. 16, namely we 
introduce a 3Nx3N tensor Vs, which projects on either 
the intramolecular modes or the ensemble of K and inter- 
molecular modes. Vs acts on the 3D eigenphonons e q „, 
such that one can define the phonon-projected matrix el- 
ements as g s = E s CPse q! ,) s • d* TO (k, k + q)/ x /2M s w qi/ . 
The resulting phonon-projected A is then 



\S,S' _ 
\v — 



1 



JV(0) N k 



9s 9*s> 



qv ■ \-/ - — k,n,m 
X (/kn — /k+q,m) <>(ek+q,r, 



£kr; 



W qi ,).(16) 



The total A is Zs,s> ^ = Es.s' E q „ V 



s,s' 



The 



contribution of each subspace S is computed as J2 S , X S ' S ' , 
where we add both the diagonal term and the usually 
very small off-diagonal contributions. The results of this 
analysis are reported in Tab. I and Fig. 6. 



TABLE I: Adiabatic A AD and non-adiabatic A computed via 
Eqs. 13 and 12, respectively, for selected electron-phonon cou- 
plings, corresponding to Fig. 6. We report also the phonon 
frequency logarithmic average o;i og for both the adiabatic and 
non-adiabatic formulations, "full el-phon" means that all 
terms are taken from the ab-initio calculation of the electron- 
phonon coupling, "local el-phon" means that that only local 
terms are retained in d qs , , (R), while "non-local el-phon" 
refers to the case where only off-site terms are taken in 
d^, n , (R). In the "local el-phon with intra phonons" not only 
the deformation potential but also the phonon eigenmodes are 
projected on the molecule. 



model 




A 


i > AD 


t^log 








(meV) 


(meV) 


full el-phon 


0.88 


0.73 


25 


18 


non-local el-phon 


0.65 


0.60 


17 


14 


local el-phon 


0.25 


0.17 


63 


41 


local el-phon with intra phonons 


0.20 


0.12 


93 


68 



(15) E 



In Fig. 6 we plot the Eliashberg function a 2 F{uj) 



q/y 



A q „w q! , 5(u - 



ujqy)/{2N q ), and the integral A(w) 



2 J Q doj'a 2 F(uj')/uj' , namely the frequency resolved 
electron-phonon coupling. The first row is the total A 
computed by means of Eqs. 13 and 12 by including all (lo- 
cal and non-local) electron-phonon contributions. This 
result has been already reported in Ref. 16. In the sec- 
ond row, we plot A where only non-local matrix elements 
are taken. This accounts for the 72% (80%) of the total 
A in the adiabatic (nonadiabatic) formulation. In this 
case, the main contribution to the electron-phonon cou- 
pling comes from intermolecular soft phonon modes, with 
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NONADIABATIC 
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0.0 

0.6 
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0.2 
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full el-phon coupling 




non-local el-phon coupling 




local el-phon coupling 




local el-phon coupling 

with intramolecular projected phonons 
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full el-phon coupling 



X((D) 
p 

a F(co) 



non-local el-phon coupling 



local el-phon coupling 



local el-phon coupling 

with intramolecular projected phonons 
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0.4 
-I 0.2 
0.0 



1200 



1500 
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1200 
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phonon frequency (cm" 1 ) 



FIG. 6: Eliashberg a 2 F and integrated A for various electron-phonon coupling models. All ab-initio elements are taken in the 
first row, the second row is for the non-local d^ s , n , (R) terms, the third and fourth rows are for the on-site coupling. In the 
latter also the phonon eigenmodes are projected on the intramolecular subspace. To the left we show quantities computed by 
the adiabatic approximation (Eq. 13), to the right those evaluated by Eq. 12. 



strong spectral weight at low frequencies (< 500cm _1 ). 
In the third row of Fig. 6, we plot the coupling aris- 
ing only from local matrix elements. The corresponding 
a 2 F(uj) is peaked around frequencies related to in-plane 
molecular phonon modes, which give the main contri- 
bution to the local electron-phonon coupling. In the 
lower row, not only the deformation potential matrix el- 
ements but also the phonon modes are projected on each 
molecule. The band structure (ek n ) and the phonon spec- 



trum (ujq U ) are instead unmodified with respect to the 
full solid. Therefore, the latter case is the closest esti- 
mate of the electron-phonon coupling of a single (doped) 
molecule placed in the crystal metallic environment. The 
projection of the phonon eigenmodes on the molecular 
subspace further reduces A. We reach therefore one of the 
main conclusions of this work. The purely molecular con- 
tributions (filtered in both the wave function and eigen- 
phonons) are such that the resulting A accounts for only 
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23% (17%) of the full A in the adiabatic (non-adiabatic) 
formulation. All the rest comes from non-local sources. 
Note also in Tab. I that the phonon frequency logarith- 
mic average uJi og varies substantially with respect to the 
model. The non-local contributions drastically reduce 
the cjiog value, that corresponds to a lower estimate of T c , 
which depends linearly on oj\ og , according to McMillan 29 . 



Technical details for the electron-phonon 
calculations 



In order to evaluate the electron-phonon coupling in 
the K3-picene, we first carried out phonon calculations 
in the density functional perturbation theory framework 
(DFPT) 30 onaiV q = 2x2x2 grid of phonon mo- 
menta q, as illustrated in Ref. 16. The electronic grid 
used in DFPT has a mesh of 2 x 2 x 2 points with a 
Methfesscl-Paxton smearing of 0.03 Ry. For each phonon 
mode v with momentum q we computed both the "nona- 
diabatic" electron-phonon interaction in Eq. 12, and the 
"adiabatic" one in Eq. 13. The k summation in both 
Equations has been performed by means of the Wan- 
nier interpolation technique 27 ' 31-33 in the MLWF basis. 
The convergence in the k summation was reached for a 
Nk — 60 x 60 x 60 momentum grid, with temperature 
and smearing given by T = 150 K and a = 4.3 meV, 
respectively, with T — 3a. 

The matrix elements g£ n k+qm have been computed 
by selecting the local and non-local contributions in the 
deformation potential (through the MLWF via Eq. 15) 
and the phonon modes (through the phonon projector 
V S ). 

With the smearing and k-point grid at convergence, 
we checked the accuracy of the extrapolated value of the 
Fermi level ep by comparing it to the value obtained 
by an ab-initio calculation with a large k-grid (N^ — 
10 x 10 x 10) and a small Gaussian smearing (0.002 Ry). 
The uncertainty in ep is of the order of 1 meV, giving rise 
to a change in the density of state N(0) of 1% only. Thus, 
a comparable error occurs in A from the ep position. 

Thanks to the Wannier interpolation the electronic k- 
summation in the electron-phonon coupling A q is con- 
verged for each q. The main residual error in the total 
A comes from the coarse A q = 2x2x2 momentum grid 
used in the q-summation of A q . An estimate of this er- 
ror was done in Ref. 16, by studying the fluctuations of 
A q over the q-point sampling. We found that the uncer- 
tainty on A goes from about 20% for the adiabatic values 
to less than 15% for the non-adiabatic estimates, while 
the error on the relative contributions (i.e. ratio of local 
A over non-local A) is even smaller (less than 10%). 



V. SCREENING OF LOCAL 
ELECTRON-PHONON COUPLINGS 

The result presented in the previous Section is com- 
pletely unexpected, if one follows the common wisdom 
that molecular crystals can be reliably described by 
molecular derived quantities 34,35 . For instance, molec- 
ular electron-phonon calculations have been used to 
compute A in molecular crystals 36,37 . This has been 
the case of some previously published works on the 
doped C6o 12,13,38 ~ 40 , and the newly discovered "aro- 
matic" superconductors 14 ' 36 ' 37,41 ' 42 . Here, we show that 
at least in the K3-picene, molecular only calculations are 
not reliable to predict the crystal total A. Presumably, 
this applies also to the whole series of new aromatic su- 
perconductors, where the physics should be similar 9 . 

To explain why the intramolecular electron-phonon 
coupling is so weak in the crystal, we carried out elec- 
tronic structure and phonon calculations for the isolated 
neutral molecule in the same Quantum-espresso 26 PW 
framework. The resulting electron-phonon coupling val- 
ues are in good agreement with previous molecular cal- 
culations by T. Kato 14 . 

The isolated molecule calculations allowed us to com- 
pute also the molecular deformation potential in the MO 
representation: 



is MOL I 



MOL 



sv s M c9 L 

5u s 



w 



MOL^ 



(17) 



where now d^^ fOL is q- independent. By replacing 
d^„(0) with d s m ™ OL in Eq. 14, and by taking only the 
local contributions (namely R = and (to, n) running 
on the same molecule), one can compare directly the dif- 

fiyMOL . 

ference between — | CF and °^ SCF on the resulting A. 

Indeed, we have already shown that the MOs w^ OL are 
very close to the molecular MLWF w m n (see Fig. 2), so 
that a difference in A can come only from the deforma- 
tion potential operator. Moreover, we noticed that the 
q-dependencc of d^ s „(0) is very weak, and so a direct 
comparison can be made at each crystal phonon momen- 

fiyMOL 

turn q. The molecular — is the "bare" one, while 
^ CF is screened by the crystal environment and by the 
partially occupied metallic bands of LUMO+1 character. 

For here on, our estimates of A are only based on the 
adiabatic approximation in Eq. 13. Although we have 
seen that it is less accurate than the nonadiabatic for- 
mula, we are going to use it because from the adiabatic 
formulation it is easier to make the connection to the 
molecular approximation formula 43 for A, and make the 
comparison with previous works (which mainly used the 
adiabatic approximation). The results are plotted in 
Fig. 7 and reported in Tab. II. By comparing panels (a) 
and (b), it turns out that the total coupling A with the 
screened deformation potential is about 4 times weaker, 
which implies that on average the electron-phonon ma- 
trix elements g^ n 

,k+qm are twice smaller than the "bare" 



10 



ones of the isolated neutral molecule. 

Therefore, we reach our second main conclusion of 
this work. Describing correctly the effect of the metallic 
screening provided by the crystal environment to the de- 
formation potential is critical to get the right estimate of 
the electron-phonon coupling. 

Now, let us analyze in details the effect of the metallic 
crystal environment on the dynamical matrix, and so on 
the phonons. In Eq. 13, we replace the phonon eigen- 
values ujqy and eigenvectors e q „, with the corresponding 
molecular ujff OL and eff OL , computed for the isolated 
undoped picene molecule. The results are reported in 
Figs. 7(c) and 7(d). If compared to panels (a) and (b), 
there is a global frequency softening of 50 cm" 1 for the in- 
plane phonons in the crystal induced by the doping. The 
second effect is a remodulation of the frequency depen- 
dence of the electron-phonon coupling strength. In the 
a F(u>) obtained with molecular phonons, the coupling is 
mostly peaked around 1600 cm" 1 , while it is much more 
broadly distributed in the crystal phonons. 

From this analysis we can conclude that upon doping 
the metallic environment provided by the crystal strongly 
affects both the deformation potential and the dynami- 
cal matrix. The metallic screening reduces the electron- 
phonon coupling strength, while it softens the phonon 
modes and makes their coupling to the charge broader in 
the phonon frequency. 

The a 2 F(oj) plotted in Fig. 7(d) for molecular phonons 
and molecular 



<5u s 



closely resembles the one published 
in Ref. 15, where the deformation potential and the dy- 
namical matrix have been computed for the undoped in- 
sulating picene crystal in the rigid doping approximation. 
Thus, the effect of the metallic screening from partially 
filled bands has been neglected in both the deformation 
potential and the dynamical matrix. The value of 
corresponding to the Eliashberg function of Fig. 7(d) is 
125 meV, very close to the value reported in Ref. 15 (126 
mcV). This is a further indication that erroneous results 
can be obtained for doped picene if the metallic screening 
is not included in the calculations. 



TABLE II: A 
tions plotted in Fig. 7. 
sequence of the Figure 



1D corresponding to the integrated a 2 F func- 



average w^g is also reported. 



Row order corresponds to the label 
The phonon frequency logarithmic 



model 




(mcv; 


crystal local el-phon with crystal intra phonons 
molecular el-phon with crystal intra phonons 
crystal local el-phon with molecular phonons 
molecular el-phon with molecular phonons 


0.20 
0.78 
0.15 
0.57 


93 
96 
110 
125 



(a) 



\((0) 

q 2 F((o) 



crystal local el-phon coupling 
with crystal intramolecular phonons 




crystal local el-phon coupling 
with molecular phonons 



(b) 



molecular el-phon coupling 

with crystal intramolecular phonons 




- (d) 



molecular el-phon coupling 
with molecular phonons 
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phonon frequency (crrf 1 ) 



FIG 



7: Eliashberg function a 2 F and the adiabatic A AD com- 



puted via Eq. 13 for both crystal local d (panels (a) and (c)) 
and purely molecular d MOL (panels (b) and (d)). In the for- 
mula, we used either the intramolecular projected phonons 
e q „ (panels (a) and (b)) or the purely molecular dynamical 
matrix (panels (c) and (d)). Note that the y axis scale of the 
left panels is 4 times wider than the scale on the right pan- 
els, to show that the magnitude of the local d (averaged over 
the phonon momenta q and the phonon modes v) is about 
4 times smaller than the magnitude of d MOL (averaged over 
the molecular phonon modes v). 



Technical details for the molecular 
electron-phonon calculations 



The molecular DFT calculations have been carried out 
with the PW basis set in the same supercell as the one of 
the K3-picene, where only one of the two molecules per 
crystal unit cell has been taken. We checked that the K 3 - 
picene supercell is large enough to get the same molecular 
levels as the ones of a much larger supercell, and thus the 
boundary effects are negligible. We left the atomic posi- 
tions of the molecule unchanged from the crystal, in such 
a way that the deformation potential calculated for the 
molecule could directly replace the one for the crystal in 
fkn,k+qm 01 Eqs. 12 and 13 without any particular ro- 
tation in the coordinate space. The DFT calculation of 
the molecule was performed in its neutral state, at the T 
point. The electron-phonon calculations were performed 
at q = (5> 3> 5) ( m crystal fractional coordinates), to 
avoid the effective charge contributions to the deforma- 
tion potential, that results in the Frohlich Hamiltonian 44 
and diverges for zone-center optical phonons. 



11 



VI. IMPACT OF DIMENSIONALITY ON 
NON-LOCAL ELECTRON-PHONON 
COUPLINGS 

In this section we want to go beyond the distinction be- 
tween intramolecular local and intermolecular non-local 
couplings, and analyze what are the most important 
electron-phonon interactions among the non-local con- 
tributions. We keep a "direct space" approach in la- 
beling the various terms, by exploiting the local picture 
provided by the molecular MLWFs. In other words, we 
aim at finding the minimal electron-phonon lattice model 
(where each lattice site represents the center of a picene 
molecule), which gives the closest possible description to 
the full "ab-initio" Hamiltonian. 

In practice, we select a subset of possible elements in 
the deformation potential matrix expressed in the molec- 
ular MLWF basis (Eq. 15). According to the set of Bra- 
vais vectors R and wannier function indexes (m, n), it is 
possible to restrict the coupling to be local (OD), unidi- 
mcnsional (ID) with molecular chains oriented along the 
b crystallographic axis, bidimensional (2D) with molec- 
ular layers spanning planes containing the a and b crys- 
tallographic axes, or the full "ab-initio" model without 
constraints. It is also possible to select the neighboring 
molecules based on their distance, therefore distinguish- 
ing between nearest neighbors (NN) and next-nearest 
neighbors (NNN) on a given direction. In Sec. Ill, we 
introduced the ID and 2D models for the hoppings. The 
same models apply also for the deformation potential ma- 
trix elements. 

The results are reported in Fig. 8 and Tab. III. The 
ID model gives 85% of the total A. The 2D model, where 
also the 2D NN and 2D NNN contributions are added, 
yields 90% of the total electron-phonon coupling. From 
the a 2 F(oj) in the upper-rightmost panel of Fig. 8, it is 
apparent that in the ID model the strongest coupling 
originates from out-of-plane vibrations and intermolecu- 
lar phonons, as a large contribution comes from frequen- 
cies below 300 cm -1 . 

This is the third important result of this work. One 
can model the system by few non-local electron-phonon 
couplings (Eq. 7) added to the local Holstein-like terms 
(Eq. 6). As the local terms are weak (see Sec. V), those 
few non-local couplings are responsible for more than 
60% of the total A. This opens the way toward an ef- 
ficient and reliable lattice modelization of the system, 
where more sophisticated many-body techniques can be 
used to deal with the electronic correlation and electron- 
phonon coupling together. 

It is striking that the main contribution to the coupling 
(more than 50%) comes from the phonon modulated 
hoppings in the b crystallographic direction. It means 
that the electron-phonon coupling is strongly anisotropic 
along molecular chains. Based on the molecular arrange- 
ment of the crystal and on the band structure, one would 
have instead expected a planar anisotropy, as the sys- 
tem is layered, with planes oriented in the a-b directions. 



On the contrary, there is no clear distinction in mag- 
nitude between the in-plane and the out-of-plane matrix 
elements, except for the hierarchy between the strong ID 
components and the rest. From this point of view, the 
material behaves more like an array of chains rather than 
an array of planes. The ID anisotropy is a consequence 
of the asymmetry of the picene armchair structure, as 
already shown in Sec. III. 

In order to generalize this argument to the experi- 
mental situations and other superconducting aromatic 
crystals, caution must be taken in view of the impor- 
tance of the arrangement of the molecules, and of dis- 
order. Even though the connection between molecular 
shape and electron-phonon anisotropy is intrinsic, and 
thus disorder- independent, disorder in experiments can 
affect the local geometry and change the total electron- 
phonon coupling strength. 




300 600 900 1200 1500 300 600 900 1200 1500 



phonon frequency (cm" 1 ) 

FIG. 8: Eliashberg function a F computed with the adia- 
batic formulation (Eq. 13) for various electron-phonon cou- 
pling models, based on a selection of the deformation potential 
d(R) terms. The ID model gives 85% of the total A, by keep- 
ing almost all of the dominant low-frequency contributions, 
related to the coupling with the intermolecular phonons. 



TABLE III: A computed via Eq. 13 for various electron- 
phonon coupling models, corresponding to Fig. 8. oj(^ is 
the phonon frequency logarithmic average. "0D" means that 
only local molecular couplings are retained in the deformation 
potential matrix elements, "ID" refers to the one-dimensional 
model of chains along the b crystallographic axis, "2D" is the 
model for the molecular herringbone layer spanning the a and 
b axes. 



model 


A AD 


(meV) 


0D 


0.25 


63 


ID 


0.74 


26 


2D 


0.77 


27 


full 3D 


0.88 


25 
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VII. CONCLUSIONS 

In this paper we went beyond what we have done in 
Ref. 16, where we carried out a detailed analysis of the 
electron-phonon coupling based on the phonon projection 
to intramolecular and intermolecular eigenmodes. Here, 
thanks to the maximally localized Wannier functions for- 
malism, we took into account also the electronic projec- 
tion of the deformation potential into local (Holstein-like) 
and non-local couplings, defined based on a molecular or- 
bital representation. We found that the purely molecular 
contribution (projected in both the deformation poten- 
tial and phonon eigenmodes) is very weak. It accounts 
for only 20% of the total A (in the adiabatic coupling for- 
mulation of Eq. 13). It is therefore impossible to predict 
the total coupling of K 3 -picene with isolated molecular 
calculations only. We understood this as an effect of the 
metallic screening which mainly reduces the deformation 
potential matrix elements. This turns out from a direct 
comparison of the projected local coupling in the crystal 
with the full coupling in the isolated molecule, carried 
out within the same theoretical framework. 

Moreover, we demonstrated that while 80% of the 
electron-phonon coupling in K3-picene is non-local, more 
than 60% of A comes just from two terms, i.e. the ID 
nearest neighbors (NN) and the ID next-nearest neigh- 



bors (NNN) phonon-modulatcd hoppings, which form 
ladder chains along the b crystallographic axis of the com- 
pound. 

The strong spatial ID anisotropy of the electron- 
phonon coupling is another interesting outcome of the 
present work. The intermolecular modes couple more 
strongly with the electrons along chains of molecules, ar- 
ranged in an ordered array of ladders. We related this 
to the picene molecular edge asymmetry in its armchair 
structure. 

Finally, we showed that the model comprising of lo- 
cal, ID NN and ID NNN contributions yields 85% of the 
total electron-phonon coupling of the crystal. This ab 
initio modelization opens the way to reliable and quanti- 
tative many-body calculations on the lattice, in order to 
study the interplay between strong electronic correlation 
(typical of a molecular crystal with fiat bands 45-47 ), and 
electron-phonon coupling in K 3 -picene. 
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